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ABSTRACT 

Supersonic  flow  past  harmonically  vibrating  two- 
dimensional  panels  and  cylindrical  shells  is  analyzed 
using  a  linearized  method-of-characteristics  procedure 
A  detailed  description  of  this  method  to  solve  the 
linearized  unsteady  potential  equation  is  given  and 
numerical  results  are  presented  to  indicate  the  nature 
of  the  aerodynamic  pressure  distributions.   Also,  com- 
parisons of  the  present  results  are  made  with  earlier 
work  by  Nelson  and  Cunningham  and  by  Anderson  which  is 
based  upon  quite  different  approaches. 
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I.   INTRODUCTION 

The  determination  of  the  pressure  distribution  over 
vibrating  panels  and  shells  is  a  prerequisite  for  the  study 
of  aeroelastic  stability.   This  problem  has  been  previously 
considered  by  a  number  of  investigators,  who  applied  either 
various  asymptotic  approximation  techniques  (e.g.,  piston 
theory)  or  operational  methods  to  the  solution  of  the  lin- 
earized unsteady  potential  equation.   For  a  recent  compre- 
hensive review  of  this  subject  refer  to  Reference  7.   In 
the  present  study  a  quite  different  approach  is  developed; 
i.e.,  the  method  of  characteristics  is  used  to  obtain  a 
solution  to  the  problem  of  supersonic  flow  past  vibrating 
panels  and  shells. 

The  major  assumptions  in  the  mathematical  model  are 
that  we  are  dealing  with  a  perfect  gas,  that  the  flow  is 
irrotational  and  supersonic,  and  that  linearization  is 
permissible.   Reference  1  was  used  to  help  develop  the 
linearized  unsteady  potential  equation  and  the  boundary 
conditions  for  panels  and  shells.   These  equations  were 
then  rewritten  in  cylindrical  coordinates. 

Steady  supersonic  flow  past  bodies  of  revolution  was 
studied  by  several  authors  using  linearized  characteristics 
methods;  i.e.,  Haack  [Ref.  9j  ,  Sauer  [kef.  lo]  ,  Oswatitsch 
and  Erdmann  [Ref.  11)  .   Teipel  [Ref .  12J  developed  a  char- 
acteristics procedure  for  two-dimensional  supersonic  flow 
past  airfoils  oscillating  at  arbitrary  frequency  using  the 
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velocity  perturbations  and  the  velocity  of  sound  perturbation 
as  dependent  variables.   Platzer  and  Sherer  [Ref.  13/  gave  an 
extension  of  the  Oswatitsch-Erdmann  method  to  slowly  oscil- 
lating bodies  of  revolution. 

In  the  present  report  a  linearized  characteristics  method 
is  presented  to  study  supersonic  flow  past  cylindrical  shells 
vibrating  at  arbitrary  frequency.   The  velocity  potential 
function  is  introduced  as  the  dependent  variable  and  the 
linearized  unsteady  potential  equation  is  solved  subject  to 
the  proper  boundary  conditions.   Although  two-dimensional 
flow  past  vibrating  panels  is  contained  in  this  solution  as 
a  special  case  for  infinite  shell  radius,  two  separate 
programs  were  written  in  Fortran  IV  code  for  the  IBM  360 
system,  one  for  panels  and  one  for  shells.   Thus,  the  panel 
program  is  complementary  to  Teipel's  work  in  that  it  uses 
the  velocity  potential  rather  than  the  velocity  perturbations 
as  dependent  variables.   The  method  of  singularities  [Ref.  4/ 
served  as  an  independent  method  to  verify  the  results. 
Anderson's  work,  on  the  other  hand,  was  used  to  compare  the 
method-of-characteristics  predictions  for  cylindrical  shells. 


II.   PROBLEM  FORMULATION 

In  nearly  all  engineering  studies  some  basic  assumptions 
have  to  be  made  to  formulate  the  problem.   In  this  study 
this  assumption   is  that  the  compressible  medium  with  which 
we   are  dealing   is  a  perfect  gas.   We  therefore  have  the 
equation  of  state  for  a  gas 

P  =  pRT  (2.1) 

Flow  of  a  perfect  gas   is  then  completely  described  by 
expressing  the  ambient  quantities  P,p  and  T  as  well  as  the 
velocity  components  u,  v,  w  as  a  function  of  position  x,  y, 
z  and  time  t.   These  six  dependent  variables  are  always 
related  by  the  conditions  of  conservation  of  mass,  momentum 
and  energy  as  developed  in  Johns   Ref .  ll . 
Conservation  of  Mass 


3 (pu)    3 (pv)    3 (pw)    3£ 
3x       3y        3z   +  3t    u  *Z,ZJ 

Conservation  of  Momentum 

1   3P    3u    -  3u  ,  Tr  3u  ,  T73u 
_  —   =  +  n  +  V +  W 

p   3x    3t      3x     3y     3z 

1   3P    3v    _  3v     3v     3v 
"p   3y-  =  It  +  U  15?  +V"3^  +W3^  (2*3) 

1   3P    3w    —  3w     3w     3w 
p   3  z    3 1      3x     3y     3  z 


/ 
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Conservation  of  Energy 


It  jj3 


+  W2 


|^(up)  +l^(vp)  +  fl(wP)]  (2.4) 


where 

W2  =  u2  +  v2  +  w2.  (2.5) 

Equation  2.4  can  also  be  stated  as  (See  Ref .  1) 


§|=0  (2.6) 

The  condition  of  conservation  of  energy  therefore  reduces 
to  the  simple  statement  that  the  entropy  of  a  fluid  particle 
remains  constant.   For  later  convenience  we  now  rewrite  the 
continuity  equation,  Eq.  2.2 


Dp     (3u    3v    3w\ 
DT  +  p  (ax  +  3^  +  3^ j 

Introducing  the  speed  of  sound 

••-if  (2.8, 

we  have 


Dp   _  3p  DP  =  1   DP 

Dt     3P  Dt  "'   2  Dt 
a 

Therefore  the  continuity  equation  can  also  be  written 


DP  =  _pa2  du         3v    3wN. 
Dt  Ux  +  3^  +Ti-) 


The  velocity  potential  is  introduced  to  reduce  the 
number  of  dependent  variables.   Its  existence  depends  upon 
the  condition  of  irrotationality  of  the  flow.   This  condition 
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is   expressed  mathematically  by  the  vanishing  of  the  curl 
of  the  velocity  vector  and  can  be  written  in  component  form 


9v    9u       9w    9v        9u    9w 

9x  -  9y  =  °'  9y  "  dz   =  °/   Tz  ~  9x  =  ° •  (2.10) 


Ref.  2 


that  flows 


It  can  be  shown  by  Kelvin's  theorem 
starting  either  from  reservoirs  or  from  parallel  stream- 
lines and  being  irrotational  initially  will  remain  so.   Only 
strong  shock  or  intense  heat  will  require  a  reformulation 
of  this  law.   However,  such  effects  are  usually  quite 
unimportant  in  aeroelastic  considerations.   Therefore,  Eq. 
2.10   is  necessary  and  sufficient  to  derive  the  velocity 
components  u,  v,  w  from  the  velocity  potential  <f>  such  that 

—     96        96        96  /-i-i-ix 

u=   9^  '  V  =  9^  '  W  =  9"f  '  (2'11) 

Therefore  the  continuity  equation  can  be  written 


p  Dt    9x    9y    9z    3x2   9y2    3z2 

In  vector  form  the  Eulerian  equations  of  motion  can  also 
be  written  as 


9v    _        _     i 

jr  +    (v  •  grad)v  =   -  ±   grad  P.  (2.13) 

P 

This  is  equivalent  in  irrotational  barotropic  flow  to 


9  w2       1 

y£  grad  <j>  +  grad   l-j- )  =  -  -  grad  P  (2.14) 
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or 


grad 


"3<J>        W2  [dp"| 

Tt+2-+       J-J=    ° 


Integration  therefore  gives 


8cf)    W2 


dP 
P 


C(t) 


(2.15) 


Moreover,  C(t)  =  ±   U2   whenever  the  remote  fluid  motion 

2 
consists  of  parallel  streamlines  with  velocity  U.   The 

partial  derivative  of  Eq.  2.15  with  respect  to  time  is 


8_ 
8t 


9<t>   w; 

at  +  2 


_3 

at 


r 


j 


dP 
P 


1  3P 

p  3t 


(2.16) 


In  rewriting  the  continuity  equation 


a2p 


[■ 


ap  ,  A  ap   .  ap   ap 


=  0 


(2.17) 


one  obtains  therefore  as  a  result  of  Eqs.  2.3  and  2.16 


I       * 


1   - 


2         TXX 


*  + 


* 


L-r^]    *yy  "l1"^'  ♦» 


2<f>x4)y 


2*    <t> 

Yxyy 


2^y<f»Z 


xy 


'xz 


*yz   + 


(2.18) 


2<J\ 


2<f>, 


2<J) 


2     *Xt 


2    yyt 


-|     »zt  -  I_«tt  =   0 


Let  a  denote  the  speed  of  sound.   For  infinitely  large 
speeds  of  sound  (incompressible  flow)  the  above  equation 
reduces  to  Laplace's  equation.   Equation  2.18  is  the 
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governing  differential  equation  for  general  nonstationary , 
nonviscous,  potential  flow  and  is  not  limited  to  small 
disturbances.   However,  because  of  its  nonlinear  nature, 
solutions  were  found  only  in  very  few  special  cases.   There- 
fore, it  is  necessary  to  introduce  the  small  disturbance 
concept.   Regarding  all  velocity  disturbances  small  in 
comparison  with  U,  a,  and  (U  -  a)  and  pressure  differences 
and  density  changes  small  in  comparison  with  main  stream 
pressure  and  density,  a  perturbation  potential  <£*  is  introduced 
so  that 

$  =  Ux  +  <p  (2.19) 

If  a  is  expanded  around  its  value  for  the  undisturbed  stream 
the  procedure  of  retaining  only  first  order  terms  leads  to 
the  following  equation  for  the  disturbance  velocity  potential 

<1-M*)fxx  +4>yy  +<>  „  -  2Sf  xt  -  i_<ptt  =  0  (2.20) 

where  c  is  the  constant  value  of  the  speed  of  sound  in  the 
uniform  stream  and  M  the  constant  free-stream  Mach  number. 
It  is  this  linearized  unsteady  potential  Eq.  2.20  which  is 
being  used  in  most  aeroelastic  analyses  as  the  basic  equation 
to  predict  the  oscillatory  pressures  and,  generalized  aero- 
dynamic forces.   However,  it  must  be  emphasized  that  its 
validity  is  restricted  to  purely  subsonic  or  supersonic 
flows,  as  well  as  to  sufficiently  unsteady  transonic  flows. 
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Expressing  Eq.  2.20  in  cylindrical  coordinates  one  has 
^xx  +  cPrr  +  -fr  +  —  fGG  ~   ±-     [cptt  +  2UflKt  +  U2^. 


r2 


xx 


=  0 


(2.21) 

There  are  two  main  boundary  conditions  applicable  to  a 
vibrating  cylinder  in  supersonic  flow.   The  first  is  that 
the  panel  surface  is  impenetrable  to  the  gas,  i.e.,  the  flow 
is  tangential  to  the  surface  at  each  instant  of  time.   The 
second  condition  is  prescribed  by  the  Sommerfeld  radiation 
condition  which  requires  waves  to  propagate  away  from  sources 
of  disturbance  toward  infinity.   Furthermore,  uniform  super- 
sonic flow  must  exist  upstream  of  the  panel  leading  edge, 
since  no  disturbances  can  propagate  upstream. 

If  the  equation  of  the  surface  of  a  body  which  is  moving 
in  a  time-dependent  manner  is  given  by 

F(x,y,z,t)  =  0  (2.22) 

then  it  can  be  shown  that  the  flow  tangency  condition 
requires  [kef.  8] 

DF  9F  _      3F  3F  9F 

Dt=3t      +u      3x"+v3y~+w3¥=0  (2.23) 

This  condition  states  that  the  rate  of  change  of  the  numerical 
value  of  the  function  F  is  zero,  when  we  follow  the  motion 
of  a  particular  fluid  element  (substantial  derivative) ,  so 
that  the  element  continually  touches  the  surface  F  =  0. 
Consider  a  two-dimensional  panel  (Fig.  la)  whose  chord  is 
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aligned  with  the  x-axis  and  whose  upper  surface  is  exposed 
to  a  uniform  flow  of  velocity  U  in  the  positive  x-direction. 
We  then  find  that  we  can  describe  the  body  by 

F  =  z-h(x,t)  =  0  (2.24) 

Since   _  =i    we  have  from  the  boundary  condition 
8z 

w  =  |^  +  u  |^      For  z=  h(x,t)  (2.25) 

Imposing  the  linear  perturbation  assumption  Eq.  2.19  we  have 

4  =  Ux  +  <p  (2.19) 

where  the  disturbance  velocity  components  are  obtained  from 

u  =  U  '        w  =  H  (2-26> 

dX  d  Z 

and  satisfy  the  order-of -magnitude  requirement 

u,  w   <<  U  .  (2.27) 

The  terms  u  — —  can  then  be  neglected  compared  to  the  much 

dX 

larger  term   U  --  .   This  leads  to 

d  X 

dh  ah 

W=it+U3^  (2-28> 

It  can  be  shown  that  it  is  possible  to  express  the  flow 
tangency  condition  in  the  plane  z  =  0.   Taylor  series 
expansion  gives 

/     *.%      t       n+    4.x    v,  9w(x,0+,t) 
w(x,z,t)  =  w(z,0  ,t)  +  h  ^ 


+  h2  32w(x,0+,t)  + + 


2«      ,_2 


(2.29) 

3zr 
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FIGURE  la   TWO-DIKENSIONAL  PANEL 
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FIGURE  lb   PROFILE  AND  CROSS  SECTION 
OF  CYLINDRICAL  SHELL 
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Using  the  same  argument  as  above  the  product  terms  are 
neglected  so  that  the  flow  tangency  condition  equals 


w  =  *h   +  U  ^ 
3t      3x 


for  z  =  0 


(2.30) 


A  similar  consideration  for  the  cylindrical  shell  (Fig.  lb) 
whose  vibration  mode  is  given  by 


h(x,6,t 


)  =  Z(x)  Cos  nQ   e 


iwt 


(2.31) 


results  in  the  following  boundary  conditions 


3^ 


r=R 


=      U 


=    0 


az 

3* 


+    iwZ 


Cos    riQ    e 


iwt 


for      0    <    x    <    L 
for      x    <    0 


(2.32) 


The  calculation  of  perturbation  pressures  and  pressure 
coefficients  begins  with  Eq.  2.15 


11  +  wf.  + 


dP 


=   C  (t) 


3t    2      |p 
For  uniform  parallel  flow  far  upstream  we  have 

Hi 

C(t)  =  2   . 

Using  again  the  perturbation  assumption,  Eq.  2.19 

4  =  ux  +<p 
and  noting  that 

2   /   v  2  i      *  2 


W2  = 


o  +  HE   +  if-  +/lf 

3X       3y      3Z 


(2.5) 


we    can   approximate 


3<j>    +      W; 
3t  2 


3t         2      +  U  3x 
18 


(2.33) 


by  dropping  second  order'  terms.   Also  the  density-pressure 
relation  is  expanded  as 


p(p)    =   p      +    dP 
dP 


p-p        +   — 


/d2p 


2!     ^P2 


P-P      + 


(2.34) 


dp 


1    + 


dp 
dp 


P-Pc 


-  fh-  —  1 

p . J  L  PY    J 


dP    4       p 


(2.35) 

Hence,  the  perturbation  pressure  is  related  to  the  perturbation 
potential  by 


p-p   =  -  p 


3        3 
It  +  U  3x" 


(2.36) 


Writing  this  equation  in  pressure-coefficient  form  one  has 


Cp  = 


P-Pco 

2 

3£ 

-    2    3j£ 

1           TT2 

U2 

3t 

U    3X 

2    P=oU 

(2.37) 


Consistent  with  the  assumed  cylinder  deflection,  Eq.  2.31  a 
velocity  potential  is  now  assumed  in  the  form 


<?  =  <Mx,r)  Cos  ne  e1(4)t 


(2.38) 


Substituting  this  into  the  Eq.  2.21  results  in 

2iwU 


(1-M2)     $         +    $         +   i    $       -    — 
vx   n    ;     *xx   t    *rr        r      r         r 


~    $x    T    a~? 


■  Si#  =  0 


(2.39) 


This  is  the  basic  equation  for  $(x,r),  which  needs  to  be 
solved,  subject  to  the  boundary  conditions 
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9_$ 
3r 


3Z    •  „ 
U  j^  +    lwZ 


r=R 


=  0 


For  0  <  x  <  L 


For  x  <  0 


(2.40) 


x      r 

When  non-dimensional  coordinates  x  =  j-f  r  =  f"  and  the 

reduced  frequency  K  =  y-  are  introduced  then  it  is  found 


1  2 

(1-M2)$ +  $    +  ~  $   -  EL  $-  2iKM2 a   +  K2M2  $  =  0 

xx    rr      ^-    r  2  ^ 


(2.41) 


The  bar  has  been  omitted  for  simplicity.   The  boundary 
condition  is 


$. 


r=R 


fe+  iKZ 


0  <_  x  <_   L 
x  <  0 


(2.42) 


=  0 


Converting  Cp  to  cylindrical  coordinates  and  using  the  above 
non-dimensional  coordinates  one  has 


Cp  =  -2 


iK$  - 


Cos  nOe  lKt 


(2.43) 
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III.   NUMERICAL  TECHNIQUE 

A.   METHOD  OF  CHARACTERISTICS 

Consider  the  following  pair  of  simultaneous  first  order 
partial  differential  equations  for  the  dependent  variables 
u  and  v, 


A    **  gu  9V 

1  ax    l  ay    l  ax 


+  Di  IX  =  o 
1  dY 


(3.1) 


a2  |H  +  b0  |H 

2  3x     2  9y 


+  C 


'v 


2  ^+  D 


1Z  =  0 


2  ay 


where  the  coefficients  A, ,  .  .  .  ,  D~  are  functions  of  u 
and  v  but  not  of  x  and  y. 


a-A/m/l/-   Coa/o/t/oa/S 


FIGURE    2 
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Assume  that  in  Figure  1  the  solution  is  known  up  to  the 
curve  CPC.   At  P  one  knows  continuously  dif f erentiable 
values  of  u  and  v  along  the  curve  CPC  as  well  as  all  their 
derivatives  in  the  directions  that  lie  below  the  curve. 
However,  the  question  needs  to  be  answered  whether  there  is 
enough  information  to  deduce  the  directional  derivatives 
of  u  and  v  at  P  in  directions  that  lie  above  the  curve. 

The  directional  derivative  of  u  in  the  direction  s  is: 


3u 

as 


3u_  3_x    3u  3y 
3x  3s  +  3y  3s 


(3.2) 

so  that  the  directional  derivative  in  any  direction  is  known 
as  soon  as  the  derivatives  3u/3x  and  3u/3y  are  known. 
To  find  the  values  of  3u/3x  and  3u/3y,  write  out 
equations  3.1  at  point  P  and  the  directional  differentials 
of  u  and  v  taken  along  CPC  at  P.   In  matrix  form  the  four 
equations  are: 


Al 

Bl 

ci 

DT 

Tu/3x 

D  " 

A2 

B2 

C2 

D2 

3u/3y 

0 

dx 

.dy 

0 

0 

3v/3x 

= 

du 

0 

0 

dx 

d^ 

3v/3y 

dv 

(3.3) 


With  u  and  v  known  at  P  the  coefficients  A-^ ,  .  .  .  .  ,  D^ 
will  also  be  known.   When  the' directions  of  CPC  are  known 
then  dx  and  dy  can  be  found.   If  u  and  v  are  then  known  along 


In  Eq.  3.3  du  is  an  abbreviation  for  3u/3s  ds,  where 
8u/9s  is  given  by  Figure  3.2  and  s  now  measures  distance  p 
along  the  curve  CPC.   The  symbols  dv,  dx,  and  dy  stand  for 
similar  abbreviations. 
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CPC,  du  and  dv  will  be  known.   Equations  3.3  therefore 
constitute  a  pair  of  four  simultaneous  linear  algebraic 
equations  for  the  four  first  derivatives.   If  the  determinant 
of  this  matrix  is  not  zero,  there  is  a  unique  solution.   In 
this  case  satisfaction  of  these  governing  equations  requires 
that  the  directional  derivatives  have  the  same  value  above 
and  below  the  line  CPC.   If  the  matrix  determinant  had 
equaled  zero,  there  would  not  be  a  unique  solution  and 
there  could  be  discontinuities  across  the  line  CPC. 

Expanding  the  determinant  of  Eqs.  3.3  and  setting  it 
equal  to  zero  one  has 

(A  C2-A2C1)dy2  -  (A  D2-A2D1+B1C2-B2C1)dxdy  + 

(B-^-B^  )dx2  =  0  (3.4) 

This  is  a  quadratic  equation  for  the  slope  dy/dx.   If  the 
direction  of  the  curve  CPC  at  P  is  such  that  it  has  a  slope 
satisfying  Eq.  3.4,  then  the  derivatives  of  u  and  v  are  not 
uniquely  determined  by  the  values  of  u  and  v  along  the  curve. 
Such  a  direction  is  called  a  characteristic  direction.   The 
quadratic  Eq.  3.4  gives  two  real  slopes,  one  real  slope  or 
a  pair  of  complex  slopes,  depending  upon  whether  the 
discriminant 

(A1D2-A2D1+B1C2-B2C1)  2  -  4  (A -.C^A^)  (B-^-B^)        (3.5) 

is  positive,  zero  or  negative.   This  is  also  the  criterion 
for  cataloging  the  system  Eq.  3.1  as  hyperbolic,  parabolic, 
or  elliptic.   The  system  is  hyperbolic  at  a  point  if  there 
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are  two  real  characteristic  directions;  it  is  parabolic  if 
there  is  only  a  single  characteristic  direction;  or  it  is 
elliptic  if  there  are  no  real  characteristic  directions  at 
the  point. 

The  foregoing  analysis  can  be  repeated  for  the  single- 
order  quasi-linear  equation, 

d2i>       ,  d2i>  d2\i> 

a  3^  +  b3xTy+  C  3Y7  "  f  (3'6> 

in  which  a,  b,  c,  and  f  are  functions  of  x,  y,  ty ,    ■*■—  t    and 
3 ip/3y -   The  characteristic  directions  are  determined  from 
the  quadratic 

ady2  -  bdxdy  +  cdx2  =0  (3.7) 

and  hence  Eq.  3.7  is  hyperbolic,  parabolic  or  elliptic  at 
a  point  according  to  whether  the  discriminant  b2-  4ac  is 
positive,  zero,  or  negative. 

To  determine  the  characteristic  curves  we  return  to  the 
pair  of  first-order  Eqs.  3.1  and  suppose  that  the  system  is 
hyperbolic  throughout  the  domain  of  interest.   At  every  point 
there  are  two  roots,  (dy/dx) a   and  (dy/dx) 3  ,  to  the  quadratic 
Eqs.  3.4.   A  curve,  which  at  each  of  its  points  has  the  slope, 
(dy/dx) a  ,  is  said  to  be  a  a  characteristic;  A  curve,  whose 
slope  is  everywhere,  (dx/dy)$,  is  said  to  be  a  3  characteristic 
Thus  there  are  two  families  of  characteristic  curves  filling 
the  domain  as  shown  in  Figure  3. 
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FIGURE    3 

If  we  are  considering  a  characteristic  direction  so  that 
the  determinant  in  Eq.  3.3  is  zero,  the  right-hand  column 
must  be  compatible  with  this  if  there  are  to  be  any  solutions 
at  all  for  the  first  derivatives;  i.e.,  when  the  right-hand 
column  is  substituted  for  any  of  the  columns  on  the  left, 
the  resulting  determinant  must  also  vanish.   Replacing  the 
fourth  column  on  the  left  with  the  column  on  the  right  and 
setting  the  determinant  equal  to  zero  leads  to 


(A1B2-A2B1)du  +   (A1C2-A2C1)  ^  -  (B^-B^) 


dv  =  0.  (3.8) 


Insert  (dy/dx)a  into  Eq.  3.8;  'it  then  becomes  an  ordinary 
differential  equation  of  u  and  v  along  the  characteristic. 
A  similar  equation  can  be  obtained  along  a  8  characteristic. 
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A  method  of  attack  has  been  outlined  for  solving  hyper- 
bolic systems:   first,  locate  the  characteristic  curves; 
second,  integrate  the  ordinary  differential  Eq.  3.8  along 
the  characteristic.   This  procedure  is  called  the  method 
of  characteristics. 

B.   FORMULATION  OF  THE  EQUATIONS  PROGRAMMED 

In  Section  II,  D,  Eq.  2.41  was  developed  as  the  partial 
differential  equation  for  $(x,r)  that  needed  to  be  solved 
subject  to  the  boundary  condition  Eq.  2.42.   If  M  >  1,  the 
preceding  equation  is  of  the  hyperbolic  type  and  possesses 
real  characteristics,  which  satisfy  the  ordinary  differential 
equation 

(M2  -  l)dr2  -  dx2  =  0  (3.9) 

Introduce  along  the  characteristics  the  arc-length  ds2=M2dr2 

then  one  has  X'  =  **   =  ^EI,  r '  =  |§  -   £  (3.10) 

ds      M  ds      M 

Note:   The  angle  between  the  characteristics  and  the  x-axis 

is   Sin  a  =  -  (3.11) 

M 

An  arbitrary  function  F(x,r)  has  the  following  derivative 
along  a  characteristic 

S-P   X'  +  F   -r'=^P^F  ±±F  (3.12) 

ds    x       r  M    x   m  r 

Denoting  the  derivative  along  the  left-running  characteristic 
with  F. ,  the  derivative  along  the  right-running  characteristic 
with  F2  one  has 
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Fl  "   M  Fx  +  m  Fr'   F2  -  "IT-  Fx  "  |  Fr 


(3.13) 


Solving  for  F  .  F   results  in 
3      x'   r 


Fr  =  f  <F1  ~  F2>'    F   = 


X    2/M2-l 


ii-   (F,  +  F2) 


(3.14) 


Call  s   the  arc-length  of  the  left-running  characteristic, 
then  one  has  for  the  second  derivatives 


dF 


dF- 


12 


Hence 


12 


ds2   ' 

21         ~   dSl 

<Fl>x 

/Mz-1          ,„    % 
'         M        "    (Fl}r 

1 
M 

(3.15) 


and 


-F 


12 


M' 


(1-M2)F    +  F^ 
xx    rr 


=  -F 
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Therefore  the  basic  equation  for  $(x,r)  was  written 


-M2S>         +   M 


12    ■    Yr     ^   "'2 


ikMj 
/M2-l 


•1i-*2 


+     |K2M2-   Jl 


or 


=    0 


(3.16) 


$ 


12     - 


*l-*2         iKM 


2rM      "v^^l 


'l    +    $2 


K2    -   pr^ 


This  equation  was  then  written  in  finite  difference  form. 
Using  the  backward  difference  form  of  the  finite  dif- 
ference method  as  developed  in  Wang   Ref.  3   one  obtained 

8x1  j  ~ 


=r  [  ~  1  -\-i   +  1  -\)      Applying  this  technique  to 
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Eq.  3.16  where  the  point's  A  and  B  were  on  the  right-running 
characteristic  (See  Fig.  3)  then 


^  (B)  -  ^  (A)    *T-*2    i 


KM 


AS 


2rM     /M^T 


I  *1+  *2  1  +  K2- 


n' 


r2M2 


$    (3.17) 


FIGURE  4 

The  values  of  $ ,  $,  and  $_  in  Eq.  3.17  were  not  specified  at 
a  point.   When  this  equation  was  first  written  these  values 
were  taken  at  A,    the  previous  point.   It  was  found  in 
comparing  the  results  in  Chapter  V  that  there  was  a  consider- 
able amount  of  error  present.   In  investigating  the  source 
of  these  errors  it  was  found  that  if  the  values  for  $,  $,  and 
$2   were  averaged  between  points  A  and  B  that  these  errors 
disappeared.   The  value  for  $  at  point  B  was  not  yet  known 
but  it  was  found  by  integrating  along  the  characteristic. 
Using  trapezoidal  rule 


$(B)  =  $(A)  +  M   $2 (A)  +  $2(B)  AS 


(3.18) 


Substituting  this  into  Eq.  3.17  where  $ (B)  arises,  and  using 
the  above  averaging  technique  for  $,,  $2  and  *  results  in 
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$X(B) 


1  - 


AS 


iKM 


4r (B)M  +  2/M2-l  AS 


$n  (A) 


1+ 


AS 


iKM 


4r(A)M  "2/M^-l 


AS 


AS 


IKM 


AS2  (2 


*2(B)   4FTBTM  +  2/M^-l  AS    T"   lK      rTT 


n2    VI  - 

B)M2JJ 


[: 


AS  iKM  AS2     /  n2 

$2(A}    L4r(A>M  +    2/M^-l    As    "    4~       \K2    "   r2(A)MT-]|     +  (3.19) 


*<A>  [K2  -Fiw]f 


+  $(A) 


K"  -   rr? 


r*l 


n    1 

B)M'J 


AS 
2 


Similarly  to  solve  for  $2  at  point  B,  the  above  procedure 
was  used  except  now  move  along  the  left  running  character- 
istic and  integrate  <£>,  along  the  left  running  characteristic 
to  find  <J>(B)  . 


*2(B) 


1  -   AS     +  iK,M     Aq 
4r(B)M    2/M2-l  Ab 


=  $o(C) 


AS 


iKM 


1+  4r(C)M"  2/M^-l  AS_ 


$X(B) 


AS      iKM        _  AS2   f  2   V 

4F(B)M    W=l  AS    J"   (K2  "  rMg)M.j 


^(c) 


AS 


iKM   _       AS2 

+   -  JT.o T   AS   -  —J— 


L4r(C)M  *  2^F-1 


Kz  -7 


n< 


r"  (C)M' 


+ 


(3.20) 


*(C) 


K/  - 


n- 


rz  (C)M: 


AS 


+  <d(C) 


K: 


n' 


r'  (B)MZ 


2 


Then  there  are  two  equations  [3.19  and  3.20J  in  two 

unknowns,  $-.(B)  and  $2(B)/  which  can  be  solved  by  using 

2 

Cramer's  ruleT   The  following  abbreviations  are  defined  as 


Crandall,  S.  H.,  Engineering  Analysis,  New  York,  1956, 
p.  36. 
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A    = 


Al    = 


A2    = 


B    = 


,    _A^_^^         iKM 


1    + 


4r(B)M   +    2y^2-i    AS 
AS  j_KM 


4r(A)M  "    27^1 


AS 


1   + 


AS 


AS 


iKM 


4r(C)M      ~    2/M2-l 


AS 


iKM 


4r(B)M  2/M2-1 


AS    -      AS 


K2    -  n 

4  rz(B)Mz 


(3.21) 


C    =    $1(A)A1+    $2(A) 


AS  iKM  ,q2. 

L      liWM  ~    2/M^-l   +  4^   (k2    - 


+       ^      *(A) 


TC2     _  n 

r"(A)M^ 


r"2lAlMr 
AS    +   \    »(A)     [K.    -   p*^]* 


D    =    $2(C)A2+    $     (C) 


_    iKM 


4r (C)M         2/M2-1 


4-       AH    .     2    _ 

4       (K  Y^{ 


C)M2    J  J 


+    2    $(c) 


K2- 


r^fc 


2      1  1 

2—  AS    +  I 

C)M2J  Z 


*(A) 


Kz- 


n' 


rz  (B)M2 


AS 


In  terms  of  these  abbreviations  the  two  equations  are 


Afc^B)  +  B$2(B)  =  C 


Bfc^B)  +  A  $2(B)  =  D 
Then  using  Cramer's  rule 

^(B)  =  CA  -  BD/A2  - 


B: 


$  (B)  =  DA  -  CD/A2  -  B2 


(3.22) 


(3.23) 


$,  was  integrated  to  find  $te)along  the  right  running  char- 
acteristic and  as  a  check  $2  was  also  integrated  along  the 
left  running  characteristic  using  the  trapezoidal  rule 
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Q1 (B)    =    $(A)    +   - 

2 


$2  (A)     +    $2  (B) 


<J>2  (B)     =    <MC)     + 


i[*i 


(C)     +    $1(B) 


AS 


AS 


(3.24) 


The   difference 


was   a  measure   of    the 


S1  (B)  -  <i>2  (B) 

accuracy.  If  it  becomes  too  large,  the  grid  size  must  be 
made  smaller.  $2  (B)  and  <J>2  (B)  were  averaged  together  for 
the  best  results. 


*(B) 


1  (B)  +  $2  (B 


3 


(3.25) 


The  above  procedure  will  suffice  for  a  general  field 
point  which  is  not  subject  to  a  boundary  condition  once  the 
procedure  has  been  initiated. 

On  the  initial  Mach  line  it  is  known  that  $  =  0,  $-j_  =  0, 
A$,  =  0,  therefore  from  Eqs.  3.20  and  3.21  (See  Fig.  5). 


A$  (E)  =  A  *  (D) 
2        2  2 


*2(E)  =  $2(D)  a~ 


(3.26) 


FIGURE  5 
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At  the  initial  point  on  the  initial  Mach  line  the  panel 
boundary  condition  also  applied,  i.e., 


r=R 

also   from  Eq.    3.14 
$„   = 


_9Z 
_9x 


+    iKZ 


■]■ 


M 


remembering  that  $,  =  0  and  combining  the  two  equations 


$, 


_2 

M 


<LE  -  iKZ 

L3X 


(3.27) 


On  the  panel  from  the  above  boundary  condition  (See 
Fig.  6) 


•(L)  -  *2(L)  =  -  ± 


H  +  iKZ 

3x 


(3.28) 


and  from  Eqs.  3.19  and  3.21  $_  (L)  +  ^(L)   was  f°und  on  the 
right  running  characteristic  and  once  again  there  exists  a 
system  of  two  equations  in  two  unknowns. 


*1(L)  -  $2<L>  =  "  M 


ii  -  iKZ 

8x 


A$  (L)  +  B$2(L)  =  C 


Using  Cramer's  rule" 


#1(L)    =     [c    +    B    2    (||  +    iKz)]/ 
*2(L)    =    ["   A  S   (ft  +    HI    /  A+ 


A+B 


(3.29) 


(3.30) 


Ibid. 
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FIGURE  6 

To  write  the  program  for  a  vibrating  panel  or  the  two- 
dimensional  case,  the  radius  of  the  cylinder  r  was  set  to 
infinity.   Inspection  of  the  basic  equations  programmed 
(Eq.  3.21)  showed  that  r  was  always  in  the  denominator, 
which  causes  these  terms  to  equal  zero.   Therefore,  terms 
containing  r  were  neglected  in  the  Fortran  IV  encoding  of 
these  equations.   For  the  vibrating  shell  or  three  dimensional 
case  the  full  equations  as  developed  were  Fortran  IV  encoded. 

Also,  it  should  be  noted  that  the  problem  was  formulated 
for  an  arbitrary  wall  deflection  Z(x).   However,  in  the 
numerical  computations  only  the  case  of  a  sinsuoidal  standing 
wave  given  by 

Z(x)  =  ZQ  Sin  (mirx)  (3.31) 

was  considered.   The  axial  wave-length  1  is  hence  equal  to 
L/2m  where  m  is  the  number  of  axial  half-waves  in  a  given 
cylinder  length  L. 
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C.   TEST  CASES 

The  first  test  case  considered  a  two-dimensional  flat 
panel.   If  the  cylinder  radius  is  set  at  a  very  large  value 
in  relation  to  the  length,  then  one  has  essentially  a  flat 
panel  in  supersonic  flow. 


R  >  >  L 

When  MFREQ  is  set  at  various  values,  NFREQ  and  K  set  at  zero, 
one  then  has  supersonic  flow  over  a  wavy  wall  as  developed 
in  Johns   Ref.  1  . 


Taking  the  wavy  wall  equations  for  $,  $   and  $x  from  Johns 
Ref.  1   and  normalizing  results  in 


*  "  -      /m^-1   Sin    (x   "    v/m2"1 
>r    =    7T    Cos      it     (x    -    /fr^I      rj 


$      =      - 


/m2-i 


Cos 


-    /m2-1         r 


3 


(3.31) 
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By  using  a  large  value  of  R  in  the  three-dimensional  case, 
it  and  the  two-dimensional  case  results  were  then  compared 
to  the  wavy  wall  results.   These  results  were  found  to  agree 
to  the  fifth  significant  figure.   This  case  was  also  used 
to  remove  the  majority  of  minor  programming  errors. 

To  remove  any  other  hidden  programming  errors  it  was 
decided  to  set  each  variable  equal  to  a  finite  value  and  to 
step  through  the  first  three  characteristic  lines  by  hand. 
The  values  were  then  compared  to  those  computed  by  the 
computer.   In  this  manner  the  rest  of  the  programming  bugs 
were  found  and  removed. 
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IV.   COMPUTER  PROGRAM  ORGANIZATION  AND  OPERATION 

A.   GENERAL  ORGANIZATION 

The  numerical  methods  discussed  in  the  foregoing 
sections  have  been  coded  as  a  FORTRAN  IV  computer  program 
for  the  IBM  3  60  system. 

The  program  consists  of  a  main  control  section  and  four 
subsections  or  subroutines.   An  algorithm  showing  the 
program  information  flow  and  the  interaction  between  the 
control  section  and  the  subroutines  was  listed  in  Figure  7. 
The  control  section  reads  in  the  input  data,  computes  the 
characteristic  grid,  computes  the  initial  condition,  and 
computes  several  constants  common  to  the  subroutines.   Then 
subroutine  COMPRX  is  called  to  compute  the  coordinates 
of  the  next  grid  point  to  be  calculated,  see  Figure  8  for 
a  diagram  of  the  characteristic  grid  used.   Two  "logical  if" 
statements  then  determine  if  the  point  in  question  is 
located  on  either  the  initial  Mach  line,  the  panel,  or  at 
a  general  field  point.   Control  is  then  passed  to  that 
subroutine  to  compute  the  velocity  potential  and  its 
derivatives  along  the  characteristics  at  that  point.   The 
coefficient  of  pressure  is  computed  along  the  panel  in 
subroutine  PANPT  and  after  the  program  has  stepped  down  to 
the  last  point  on  the  last  characteristic,  a  listing  of  the 
coefficient  of  pressure  at  each  panel  point  is  printed  out. 
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FIGURE  7 
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B.   MAIN  CONTROL  PROGRAM ' OPERATION 

The  main  program  dimensions  the  storage  requirements 
in  common  block  form;  then  it  reads  three  input  data  cards. 
The  first  card  has  the  date  of  the  run  in  3A4  format.   The  .  - 
second  card  reads  in  the  fineness  of  the  grid,  number  of 
this  run,  circumferential  mode  number,  and  the  axial  mode 
number  in  415  format.   The  third  card  reads  Mach  number, 
the  reduced  frequency  number,  and  the  cylinder  radius  in 
3F10.4  format.   Input  is  then  written  to  verify  input  data. 

The  main  program  then  computes  the  Mach  angle,  step 
sizes,  and  zeros  the  counters.   The  values  for  the  velocity 
potential  and  its  derivatives  and  the  coordinates  of  each 
point  are  stored  as  a  n  x  2  matrix.   As  the  program  steps 
from  characteristic  line  to  characteristic  line,  only  those 
values  for  the  line  being  computed  and  the  line  just 
computed  are  in  storage.   The  initial  values  at  the  first 
point  on  the  initial  characteristic  is  computed  using  Eq. 
3.27  and  written  out  to  aid  in  the  program  check-out.   The 
line  counter,  "Kount" ,  is  then  set  at  two  and  subroutine 
COMPXR  is  called  to  compute  the  values  of  x  and  r  at  the 
next  point.   The  main  program  next  computes  various  quantities 
dependent  on  x  and  r  which  are  needed  in  the  other  subroutines 
The  values  for  x,  r,  and  dependent  quantities  are  printed  out 
as  a  program  debug  aid.   Two  "logical  if"  statements  next 
determine  if  the  point  under  consideration  is  on  the  initial 
Mach  line,  at  a  general  field  point,  or  on  the  panel.   The 
appropriate  subroutine  is  then  called  to  solve  for  the 
velocity  potential  and  its  derivatives  along  the 
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characteristics  at  that  point.   The  next  "logical  if" 
statement  determines  if  the  program  must  switch  lines 
or  step  along  the  present  line.   Two  more  "logical  if" 
statements  are  used  to  determine  if  this  solution  is 
completed  and  if  there  is  more  data  to  be  read  and  more 
solutions  to  be  computed.   If  "numrun"  is  set  at  two, 
the  program  will  look  for  more  data;  if  "numrun"  equals 
one,  the  program  will  terminate.   Before  it  stops  or 
reads  more  data  the  coefficient  of  pressure  along  the 
panel  will  be  printed  out  for  each  panel  point  and  an  in- 
line graph,  using  the  IBM  scientific  library  program  POLTP , 
will  plot  the  real  and  imaginary  parts  of  the  coefficient 
of  pressure. 

C.   SUBROUTINE  COMPXR  OPERATION 

Information  is  passed  between  the  calling  program  and 
the  subroutines  by  means  of  the  common  storage.   The  first 
step  in  COMPXR  is  to  determine  if  the  point  to  be  computed 
is  on  the  initial  Mach  line  or  not.   The  coordinates  of  the 
last  point  computed  are  used  and  the  step  distance  in  the 
x  and  r  directions  either  added  to  or  subtracted  from  them 
to  find  the  new  x  and  r  coordinates.   The  addition  or  sub- 
traction depends  upon  whether  we  step  up  to  a  new  character- 
istic line  or  down  the  present  characteristic  line. 
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D.  SUBROUTINE  MACHLN  OPERATION 

In  subroutine  MACHLN  the  boundary  conditions  of  the 
equations  developed  in  the  foregoing  section  are  enforced 
such  that  the  velocity  potential  and  its  derivatives  along 
the  right  running  characteristic  are  set  to  zero.   Equations 
3.26  were  coded  in  FORTRAN  to  solve  for  $   at  the  point  in 
question.   The  value  of  $2  and  the  variables  used  to  solve 
for  $2   are  printed  out  for  the  first  12  characteristic  lines 
to  aid  in  the  program  checkout  and  then  control  is  returned 
to  the  main  program. 

E.  OPERATION  OF  SUBROUTINE  PANPT 

In  subroutine  PANPT  the  boundary  condition  along  the 
panel  is  enforced.   The  velocity  potential  and  its  derivatives 
along  the  characteristics  are  computed  at  the  point  in 
question,  using  Eq.  3.29.   The  real  and  imaginary  parts  of 
the  coefficient  of  pressure  are  stored  for  future  read  out. 
Control  is  then  returned  to  the  main  program. 

F.  OPERATION  OF  SUBROUTINE  GENFPT 

GENFPT  computes  the  velocity  potential  and  its  first 
derivatives  along  the  right  and  left  running  characteristics. 
It  draws  the  necessary  information  from  common  storage  and 
uses  the  Eqs.  3.21,  3.22,  and  3.23,  which  were  developed  in 
the  foregoing  chapter.   To  aid  in  program  check  out  the 
variables  A,  A,,  A2 ,  B,  C,  and  D  are  printed  out  for  the 
first  12  characteristic  lines.   Control  is  then  returned  to 
the  main  program. 
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V.   RESULTS  AND  DISCUSSION 

The  case  of  a  vibrating  two-dimensional  panel  was 
chosen  as  the  first  test  of  the  accuracy  of  the  character- 
istics method.   These  results  were  compared  to  the  results 
obtained  by  evaluating  the  method-of-singularities  integrals 
developed  in  Ref.  4  for  flutter  of  a  two-dimensional  panel 
with  one  surface  exposed  to  supersonic  flow.   These  integrals 
were  evaluated  by  applying  a  12-point  Gaussian  Quadrature. 
The  results  of  these  evaluations  were  available  in  unpublished 
notes  by  M.  F.  Platzer. 

The  comparisons  for  M  =  v^T~at  a  reduced  frequency  of  0.5 
and  an  axial  mode  number  of  one  are  shown  in  Figures  8-11. 
Figures  8  and  9  are  plots  of  the  real  and  imaginary  parts  of 
the  coefficient  of  pressure  obtained  by  the  method  of  singu- 
larities.  Figures  10  and  11  are  plots  of  the  real  and 
imaginary  parts  of  the  Cp  for  the  characteristics  method. 
When  these  tests  were  first  run  it  was  found  that  the  real 
parts  of  the  Cp  plots  were  identical,  as  shown  in  Figures  8 
and  10.   The  imaginary  parts  of  Cp,  shown  in  Figures  9  and  11, 
did  not  agree.   In  rechecking  the  formulation  of  the  equations 
it  was  found  that  the  values  of  $  ,  $  and   $2  nad  not  been 
averaged  as  explained  in  Chapter  II,  Section  B,  where  Eqs. 
3.19  and  3.20  were  formulated.   When  these  average  values 
were  incorporated  into  the  equations,  the  computer  program 
rewritten  and  the  program  rerun,  then  the  two  plots  were 
found  to  be  identical  for  both  the  real  and  imaginary  parts 

of  the  Cp. 
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Figures  12-14  are  plots  for  a  reduced  frequency  of  two, 
an  axial  mode  of  four,  and  the  same  Mach  number.   These 
comparisons  of  the  real  and  imaginary  parts  of  the  Cp  were 
also  found  to  be  identical. 

In  attempting  to  verify  the  results  of  the  vibrating 
shell  program,  Ref.  5  was  used.   Reference  5  includes  a 
study  of  supersonic  potential  flow  past  an  infinitely  long 
vibrating  cylinder.   It  was  not  known  what  effect  the  dif- 
ference in  an  infinitely  long  cylinder  vice  a  finite  length 
cylinder  would  have  on  the  outcome,  but  these  were  the  only 
published  results  found  that  could  be  used  for  a  comparison 
of  the  vibrating  shell  program. 

Two  parameters  were  used  to  compare  the  results  of  the 
two  studies:   the  maximum  amplitude,  A,,  of  the  sinsuoidal 
Cp  wave,  and  the  phase  angle,  ty±,    of  this  wave.   Equation  14 
in  Ref.  5  is  used  to  equate  Cp  to  the  value  of  A,  listed  in 
Figure  2a  of  Ref.  5. 

_  zo  pU2      tuttR  /mux         . 

p(x,R,S)    =rif     —  Cos   n6   AX    Cos    1-^-  +    <J>; 


p   is  defined  as  p  =  p  -  pro     where 


P  -  P°° 


(5.1) 
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and  the  terms,  Z   l,  and  Cos(x)  were  all  taken  as  equal 


to  one  then 


_  p       ? 

P     h    pU^      M     1 


M 
Al  "  Cp   2m  tt  (5-2) 


Three  different  values  of  n  were  computed  and  the  coefficients 
of  pressure  were  listed  in  Figure  17  as  a  test  case  for  this 
comparison.   For  the  case  n  =  0,  (See  Figure  2a  of  Ref.  5) 
we  find  that  for  the  ratio  of  axial  wave  length/cylinder 
radius  £/R  =  3.33,  A   equals  1.05.   Using  Eq .  5 . 2  to  determine 
A   for  a  maximum  Cp  of  5.63  in  Figure  17,  A-i  for  this  study 
would  equal  1.044.   Using  this  same  procedure  for  n  =  4  and  8 

n  =  4    A  (ref  5)  =  1.30    A   =  1.445 
n  =  8     A,  (ref  5)  =  1.80    A.,  =  1.85. 

To  compare  the  phase  angle  ip -,,  Figure  2b  in  Ref.  5  was 
entered  for  I /R   =  3.33  and  the  three  values  of  n  taken 
previously.   This  shift,  i|> -.,  was  calculated  from  results 
(Fig.  17)  by  determining  the  zero  crossing  point. 

n  =  0  59.8 
n  =  4  "  57.7 
n  =  8       48.9 

Converting  these  zero  crossing  points  to  angular  shifts, 
it  was  found  that 
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n  =    0 

$1   -  4.8 

*!    (ref    5) 

=    5 

n   =    4 

ip,    =    13.5 

*1    (ref    5) 

=    10 

n  =    8 

i|>      =    55.4 

*x    (ref    5) 

=    60 

which  compare  favorably  with  Anderson's  values  listed  in 
the  third  column. 
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COEFFICIENT  OF  PRESSURE  COMPUTER  OUTPUT  LISTING 
FOR  CASE  OF  VIBRATING  SHELL 

M=  3.5,  m=  3,  K=  0,  r  =  0.20  and  n  =  0,  1,  3, 
where  L  is  divided  into  120  grid  points. 


Grid 

Point 

n  =  0 

n  =  4 

n  =  8 

1. 

5.6198 

5.6198 

5.6198 

2. 

5.5677 

5.5538 

5.5642 

3. 

5.4813 

5.4266 

5.4676 

4. 

5.3613 

5.2396 

5.3307 

5. 

5.2083 

4.9950 

5.1546 

6. 

5.0233 

4.6954 

4.9404 

7. 

4.8075 

4.3440 

4.6897 

8. 

4.5621 

3.9444 

4.4041 

9. 

4.2887 

3.5007 

4.0857 

10. 

3.9890 

3.0172 

3.7366 

11. 

3.6647 

2.4988 

3.3591 

12. 

3.3180 

1.9504 

2.9559 

13. 

2.9509 

1.3775 

2.5295 

14. 

2.5657 

0.7855 

2.0828 

15. 

2.1648 

0.1801 

1.6188 

16. 

1.7507 

-0.4329 

1.1407 

17. 

1.3258 

-1.0477 

0.6515 

18. 

0.8929 

-1.6586 

0.1546 

19. 

0.4546 

-2.2597 

-0.3467 

20. 

0.0136 

-2.8453 

-0.8492 

21. 

-0.4273 

-3.4100 

-1.3495 

22. 

-0.8656 

-3.9435 

-1.8442 

23. 

-1.2985 

-4.4557 

-2.3301 

24. 

-1.7232 

-4.9267 

-2.8040 

25. 

-2.1373 

-5.3572 

-3.2626 

26. 

-2.5380 

-5.7431 

-3.7028 

27. 

-2.9230 

-6.0806 

-4.1218 

28. 

-3.2899 

-6.3665 

-4.5167 

29. 

-3.6365 

'  -6.5980 

-4.8848 

30. 

-3.9605 

-6.7726 

-5.2237 

31. 

-4.2600 

-6.8886 

-5.5309 

32. 

-4.5331 

-6.9446 

-5.8044 

33. 

-4.7782 

-6.9398 

-6.0423 

34. 

-4.9938 

-6.8738 

-6.2428 

35. 

-5.1785 

-6.7468 

-6.4046 

36. 

-5.3312 

-6.5595 

-6.5263 

37. 

-5.4509 

-6.3132 

-6.6071 

38. 

-5.5369 

-6.0095 

-6.6462 

39. 

-5.5887 

-5.6507 

-6.6432 
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Grid 

Point 

n  =  0 

n  =  4 

n  =  8 

40. 

-5.6059 

-6.5979 

-5.2394 

41. 

-5.5885 

-6.5104 

-4.7787 

42. 

-5.5366 

-6.3810 

-4.2721 

43. 

-5.4504 

-6.2105 

-3.7235 

44. 

-5.3305 

-5.9996 

-3.1371 

45. 

-5.1777 

-5.7495 

-2.5175 

46. 

-4.9929 

-5.4616 

-1.8695 

47. 

-4.7773 

-5.1375 

-1.1981 

48. 

-4.5320 

-4.7791 

-0.5086 

49. 

-4.2588 

-4.3884 

-0.1934 

50. 

-3.9592 

-3.9678 

-0.9027 

51. 

-3.6351 

-3.5196 

-1.6137 

52. 

-3.2886 

-3.0466 

-2.3205 

53. 

-2.9216 

-2.5516 

-3.0178 

54. 

-2.5366 

-2.0375 

-3.7000 

55. 

-2.1359 

-1.5074 

-4.3616 

56. 

-1.7219 

-0.9644 

-4.9975 

57. 

-1.2972 

-0.4119 

-5.6025 

58. 

-0.8645 

0.1467 

-6.1719 

59. 

-0.4263 

0.7082 

-6.7010 

60. 

0.0145 

1.2692 

-7.1856 

61. 

0.4554 

1.8261 

7.6218 

62. 

0.8935 

2.3757 

8.0060 

63. 

1.3262 

2.9145 

8.3350 

64. 

1.7508 

3.4393 

8.6062 

65. 

2.1647 

3.9467 

8.8171 

66. 

2.5653 

4.4338 

8.9659 

67. 

2.9502 

4.8975 

9.0512 

68. 

3.3170 

5.3350 

9.0720 

69. 

3.6633 

5.7435 

9.0280 

70. 

3.9872 

6.1205 

8.9191 

71. 

4.2866 

6.4636 

8.7458 

72. 

4.5596 

6.7708 

8.5092 

73. 

4.8046 

7.0400 

8.2107 

74. 

5.0200 

7.2697 

7.8522 

75. 

5.2046 

7.4582 

7.4361 

76. 

5.3571 

7.6045 

6.9652 

77. 

5.4767 

7.7075 

6.4427 

78. 

5.5626 

7.7666 

5.8721 

79. 

5.6142 

7.7813 

5.2573 

80. 

5.6313 

7.7514 

4.6027 

81. 

5.6138 

■   7.6771 

3.9126 

82. 

5.5617 

7.5588 

3.1920 

83. 

5.4755 

7.3971 

2.4459 

84. 

5.3555 

7.1929 

1.6793 

85. 

5.2025 

6.9474 

0.8978 
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Grid 

Point 

n  =  0 

n  =  4 

n  =  8 

86. 

5.0176 

6.6619 

0.1068 

87. 

4.8018 

6.3383 

-0.6882 

88. 

4.5565 

5.9783 

-1.4817 

89. 

4.2831 

5.5841 

-2.2681 

90. 

3.9834 

5.1580 

-3.0419 

91. 

3.6592 

4.7026 

-3.7977 

92. 

3.3125 

4.2205 

-4.5302 

93. 

2.9455 

3.7146 

-5.2342 

94. 

2.5604 

3.1879 

-5.9048 

95. 

2.1595 

2.6436 

-6.5374 

96. 

1.7454 

2.0848 

-7.1275 

97. 

1.3206 

1.5150 

-7.6710 

98. 

0.8878 

0.9376 

-8.1640 

99. 

0.4495 

0.3559 

-8.6031 

100. 

0.0085 

-0.2264 

-8.9853 

101. 

-0.4324 

-0.8061 

-9.3078 

102. 

-0.8706 

-1.3795 

-9.5684 

103. 

-1.3034 

-1.9434 

-9.7652 

104. 

-1.7282 

-2.4942 

-9.8968 

105. 

-2.1421 

-3.0288 

-9.9623 

106. 

-2.5429 

-3.5439 

-9.9611 

107. 

-2.9278 

-4.0364 

-9.8932 

108. 

-3.2947 

-4.5035 

-9.7591 

109. 

-3.6412 

-4.9423 

-9.5595 

110. 

-3.9652 

-5.3503 

-9.2957 

111. 

-4.2646 

-5.7250 

-8.9696 

112. 

-4.5378 

-6.0642 

-8.5833 

113. 

-4.7828 

-6.3660 

-8.1393 

114. 

-4.9984 

-6.6285 

-7.6406 

115. 

-5.1830 

-6.8502 

-7.0906 

116. 

-5.3356 

-7.0299 

-6.4928 

117. 

-5.4553 

-7.1665 

-5.8513 

118. 

-5.5413 

-7.2593 

-5^1703 

119. 

-5.5931 

-7.3078 

-4.4545 

120. 

-5.6103 

-7.3117 

-3.7084 
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VI.   SUMMARY  AND  CONCLUSIONS 

A  linearized  method  of  characteristics  procedure  was 
developed  to  compute  supersonic  flow  past  vibrating  panels 
and  shells.   Good  agreement  was  obtained  when  comparing 
with  previous  results  by  Nelson  and  Cunningham  [kef.  4j 
and  Anderson  [kef.  5_j  ,  who  used  quite  different  approaches 
to  this  problem. 

However,  more  work  needs  to  be  done  to  further  verify 
the  cylindrical  shell  computations.  In  particular,  the 
computations  of  the  generalized  aerodynamic  forces  presented 
by  Dowell  and  Widnall  [kef.  6}  should  be  verified.  For  this 
purpose  the  limiting  case  of  vanishing  shell  radius  (slender 
body  theory)  must  be  studied.  This  requires  a  reformulation 
of  the  flow  boundary  condition. 

Also,  it  should  be  noted  that  the  characteristics  method 
allows  the  incorporation  of  quite  arbitrary  vibration  modes 
(not  restricted  to  sinsuoidal  modes)  and  hence  should  prove 
to  be  a  quite  versatile  tool  for  aeroelastic  analyses  in  the 
supersonic  flow  regime. 
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APPENDIX   A:       VIBRATING    PANEL    COMPUTER   PROGRAM 
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APPENDIX  C:   VIBRATING  PANEL  SAMPLE  COMPUTER  OUTPUT 
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COEFFICENT  OF  PRESSURE 

1.0  0.0  0.0 

2.0  24.55586270  0.00408154 

3.0  22.85066452  0.01351091 

4.0  20.09254127  0.03334298 

5.0  16.40344123  0.06806011 

6.0  11.94654108  0.12131901 

7.0  6.91909287  0.19573567 

8.0  1.54377602  0.29271826 

9.0  -3.94106736  0.41235615 

10.0  -9.29188676  0.55337020 

11.0  -14.27060977  0.71312712 

12.0  -18.65509318  0.88771781 

13.0  -22.24888160  1.07209644 

14.0  -24.88984639  1.26027485 

15.0  -26.45732668  1.44556368 

16.0  -26.87745992  1.62085036 

17.0  -26.12646822  1.77890147 

18.0  -24.23175624  1.91267680 

19.0  -21.27077195  2.01564102 

20.0  -17.36767944  2.08205951 

21.0  -12.68798847  2.10726521 

22.0  -7.43137495  2.08788457 

23.0  -1.823C0542  2.02201259 

24.0  3.89625584  1.90932880 

25.0  9.48032952  1.75114890 

26.0  14.68863637  1.55040919 

27.0  19.29657652  1.31158411 

28.0  23.10531564  1.04053992 

29.0  25.95045047  0.74433023 

30.0  27.70917520  0.43094193 

31.0  28.30563578  0.10900179 

32.0  27.71423823  -0.21254394 

33.0  25.96076634  -0.52476371 

34.0  23.12126002  -0.81907450 

35.0  19.31870316  -1.08756155 

36.0  14.71766553  -1.32327050 

37.0  9.51713284  -1.52045979 

38.0  3.94183787  -1.67480302 

39.0  -1.76752874  -1.78353270 
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40.0  -7.36480060  -1.84551983 

41.0  -12.60905291  -1.86128623 

42.0  -17.27508720  -1.83294939 

43.0  -21.16322567  -1.76410272 

44.0  -24.10798772  -1.65963661 

45.0  -25.98527008  -1.52550839 

46.0  -26.71771716  -1.36847141 

47.0  -26.27804790  -1.19577514 

48.0  -24.69019390  -1.01484935 

49.0  -22.02320014  -0.83298627 

50.0  -18.41293660  -0.65703434 

51.0  -14.00676536  -0.49311690 

52.0  -9.00639664  -0.34638770 

53.0  -3.63424668  -0.22083373 

54.0  1.87132446  -0.11913365 

55.0  7.26646797  -0.04257745 

56.0  12.31253691  0.00894953 

57.0  16.78654483  0.03691793 

58.0  20.49093494  0.04404391 

59.0  23.26223271  0.03412212 

60.0  24.97820347  0.01181008 
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APPENDIX  D:   VIBRATING  SHELL  SAMPLE  COMPUTER  OUTPUT 
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